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Simple  upwind  schemes  are  the  basis  of  many  codes  of  practical 
importance  for  design  studies  in  heat  transfer  and  fluid  flow 
applications,  e.g.  the  FL0W3D  code.  [  1  With  care,  it  is  possible  to 
produce  codes  that  can  give  bounded  answers  regardless  of  raeshyspacing  or 
time-step  size.  This  note  highlights  a  rather  disturbing  inaccuracy  that 
can  arise  if  the  mesh/spacing  is  taken  too  large.  n--iv  <rJ-'  f 


The  relevant  dimensionless  measure  of  in  for  fluid  flow 

calculations  is  the  mesh  Reynolds  number  R  ,  defined  like  the  usual 

m 

Reynolds  number  R,  but  with  dx  as  length-scale  It  has  become 

customary  [2,3]  to  study  the  effects  of  large  R  using  the  viscous 

Burgers'  equation  as  a  model.  Following  the  work  of  Cheng  and  Shubin 
[2],  hereinafter  referred  to  as  "CS",  we  consider  steady  state  flows  u(x) 
which  must  satisfy 


subject  to  boundary  conditions  u(  ±1 )  =  +1  ■ 

The  boundary  conditions  of  oppositely  directed  flow  lead  at  large  R 
to  the  formation  of  a  viscous  layer  at  x  *  0  where  the  flow  u  abruptly 
changes  direction.  The  problem  is  sufficiently  simple  that  an  analytic 
solution  i3  available,  namely 

u(x)  »  -  a  tanh(  aRx/2 ) ,  (2) 

where  a  is  the  positive  root  of 


a  tanh(  <xR/2)  =  1 


-  1 


Evidently  the  viscous  layer  has  thickness  1/R;  when  this  is  small 
relative  to  the  mesh-spacing  employed  in  a  numerical  calculation,  we  shall 
see  that  troublesome  errors  may  occur. 

For  the  purposes  of  this  note,  the  antisymmetry  about  x  =  0  means 
we  need  only  consider  what  happens  on  the  interval  x  =  0  to  x  =  1 .  We 
introduce  a  uniformly  spaced  mesh,  with  points 

x^  *  i  Ax,  i  »  0,  ...,  N,  (NAx  »  1),  (4) 

and  seek  values  u^  i  »  0,  N,  to  approximate  u(x  ).  The 

operators  uAu/dx  and  d^/dx2  must  also  be  di3cretised. 

For  u3u/dx  we  use  the  customary  conservative  upwind  difference 
scheme  as  proposed  by  Patankar  [4]  for  advection  of  a  scalar,  such  as 
temperature  T.  In  that  situation  it  is  convenient  to  use  staggered 
meshes,  i.e.  u^  and  are  defined  at  points  separated  by  Ax/2.  If 

this  formalism  is  applied  to  Burgers'  equation,  it  follows  we  need  to  have 
some  way  of  approximating  u  at  points  *i+1/,2  ■  (i+1/2)Ax.  (A  slightly 
different  treatment  is  required  to  give  the  conservative  upwind  scheme 
described  in  CS.)  Inevitably  we  take 


Vl/2  “  (Ui  +  V,)/2' 


giving  (since  ui  <  0  for  i  >  0)  the  discrete  analogue  of  Burgers' 
equation 


2  Ax 


i  =■  1 


(  Ax) 
N-1. 


(6) 


(We  have  taken  the  customary  3  point  centred  difference  formula  for  the 
diffusion  operator.)  Eq.(6)  can  be  used  to  formulate  the  CS  problem  as 


i+1 


-  u l  +  ui(ui+1  -  u^,)  =  l-lu.^  -  2u.  +  u.  +  1  ), 


1, 


N  -  1, 


(7) 


subject  to  u 


0, 


-1. 


Following  CS,  we  --sum  over  equations  (7);  because  the  scheme  is 
conservative,  internal  contributions  cancel  pairwise,  giving 


”u  1 


a  a  + 
o  1 


u  u 
N-1  N 


4  r 

—  Krui 


N-1 


(8) 


Since  u^  and  u^  are  given,  this  is  a  quadratic  equation  for  u1  in 

terras  of  u  .  Writing  u  =  -  1  +  e,  we  obtain 
N-1  N-1 

u,  =-±  2-C-^  }.  c 

R  R  R 

mm  m 

By  Taylor  expanding  the  exact  solution  (2)  about  x  =  1,  we  find 


2  sinh  2(  aR/2  ) 


For  R  >>  1,  (3)  implies  a  =  1,  hence  both  e  and  e/K  are  small, 

m 

and 


Ml) 


i.a.  as  R  -*•<“,  u.  =  -  /2  (we  take  the  negative  root  because  (7) 
m  1 

assumes  <  0,  all  i  )• 

The  exact  solution  (2)  is  monotone  and  for  large  R,  u  =>  -  1 

everywhere  outside  a  narrow  layer.  Hence  the  scheme  (6)  produces  a  result 

in  error  by  =40%  if  R  is  too  large,  or  equivalently  >  1/R,  i.e. 

ra  ~ 

the  viscous  layer  is  not  resolved.  The  worrying  feature  is  that  the 

misleading  result  is  obtained  for  any  sufficiently  high  R  ,  thus  halvinq 

cn 

Ax  need  not  change  the  value  of  the  velocity  at  the  first  interior 
mesh-point.  Only  by  checking  the  spatial  distribution  of  the  computed  u. 
can  the  error  be  detected,  see  Fig.  1. 

There  remains  the  question  of  error  behaviour  at  smaller  R  .  Fig.  2 

m 

plots  the  difference  E  between  (eq.(11)  )  and  the  true  value  of 

u( Ax)  from  (2).  We  see  that  the  error  E  passes  through  zero  for 

Rm  3  3.3/  and  thus  that  E  <3%  provided  R^  <  3.8.  Over  this  range 

of  R  /  the  Patankar  scheme  outperforms  the  conservative  upwind  scheme  in 
m 

CS  which  gives  errors  of  up  to  about  17%.  Misleading  results  are  obtained 

only  at  higher  R  . 

ra 

Lest  it  be  thought  that  the  model  problem  is  too  idealised  to  be  of 

practical  use,  we  remark  that  a  40%  error  in  velocity  will  lead  to  a 

similar  error  in  heat  fluxes,  other  things  being  equal.  ft  mesh  refinement 

study  was  conducted,  using  the  FLOW3D  code,  for  the  Rayleigh-Benard 

problem  in  a  closed  box  of  aspect  ratio  3/2  (see  [5  ]  for  more  details). 

The  results  obtained,  for  a  Rayleigh  number  Ra  =  120000  and  Prandtl  number 

Pr  »  2.5  are  listed  in  Table  I.  Observe  that  the  error  in  heat  flux  Nu 

depends  on  R^  much  as  we  should  expect  from  Fig. 2,  and  that  the  error  is 

as  large  as  25%  when  R  =  7 . 

m 

Detailed  inspection  of  the  flow  field  from  the  finer  mesh  reveals  the 
presence  of  small  counter-circulating  "Moffatt"  eddies  in  the  corners  of 
the  box.  Overall,  flow  is  more  vigorous  on  the  finest  mesh,  R  =  49 


against  R  =  40  on  the  coarsest,  but  in  a  corner,  at  e.g  (x,z)  =  (17/12, 

1  '6)  we  have  that  the  vertical  component  of  velocity  U  =  -  12  as 

against  0  =  -  22  on  the  coarser  mesh  with  Ax  =  1/6.  The  greater 
z 

speeds  in  all  four  corners  are  apparently  ultimately  responsible  for  the 
increased  heat  flux  Nu,  since  the  breakdown  of  Nu  into  conductive  and 
convective  contributions  (at  z  =  1/24)  is:  4.01  and  0.96  on  this  coarse 
mesh,  3.84  and  0.12  on  the  finest  mesh.  Obviously  the  situation  is  not 
as  clear  cut  as  could  be  desired.  Studies  even  at  low  Prandtl  number  show 
that  the  linkage  between  velocity  and  temperature  fields  persists  (i.e.  we 
cannot  find  an  accurate  temperature  distribution  unless  the  velocity  is 
correspondingly  well  represented).  However,  a  two-dimensional  version  of 
the  above,  spurious  velocity  enhancement  mechanism,  must  be  a  strong 
candidate  to  explain  the  faster  flow  on  the  coarser  meshes. 
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Table  I 


Mesh  refinement  study  using  FL0W3D  for  the  Rayleigh-Benard  problem  in  a 

closed  box  of  aspect  ratio  3/2  at  a  Rayleigh  number  of  120  000  (Prandtl 

number  2.5).  R  is  mesh  Reynolds  number  and  Nu  is  the  total  heat  flux, 
m 

across  a  plane  near  a  horizontal  boundary,  normalised  with  respect  to  the 
conductive  flux  in  the  absence  of  convection. 
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